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The Go model is extended to the case when the non-native contact energies may be either attractive 
or repulsive. The folding temperature is found to increase with the energy of non-native contacts. 
The repulsive non-native contact energies may lead to folding at T = for some two-dimensional 
sequences and to reduction in complexity of disconnectivity graphs for local energy minima. 

PACS numbers: 71.28.-hd, 71.27.-Fa 



Functionally useful proteins are sequences of 
amino acids that fold rapidly under appropriate con- 
ditions (temperature range, acidity of the water solu- 
tion etc) into their native states commonly assumed 
to be their ground state configurations The dy- 
namics of folding is akin to motion in a rugged free 
energy landscape ||] and it crucially depends on two 
factors: the interactions between the amino acids 
and the target conformations. In this paper, we 
focus on the role of non-native contact energies in 
the folding process. We study this issue in a sim- 
ple model which is an extension of the so called Go 
model proposed by Go and Abe |^ 

Specifically, we consider the standard two- 
dimensional lattice model of 16 monomers. Its 
Hamiltonian is as follows 



H 



i<j 



(1) 



where Ay = 1 if monomers i and j are in contact 
and Aij = otherwise (monomers i and j are con- 
sidered to be in contact if they are seperated by one 
lattice bond and \i — j\ ^ 1). The quantity aij = —1 
if monomers i and j are in contact in the native con- 
formation and aij = a otherwise. The sequence is 
thus defined by the native conformation. We allow 
a to be attractive {a < 0) or repulsive (a > 0). In 
what follows a is assumed to be larger than -1.0 so 
that the native state is guaranteed to be a maximally 
compact conformation. In the original Go model |^ 
a = 0. It should be noted that using model (1) one 
can monitor the effect of non-native contact ener- 
gies by varying only one parameter a. Furthermore, 
the two-dimensional model is simple enough to study 
the effect of the target conformations on the folding 
dynamics. 

It should be noted that there have already been 
several studies of models with repulsion in the non- 
native contacts: a Gaussian model M, a designed 



model g], and the so called IIP+ model The 
generalized Go model that we study here allows one 
to vary the strength, and sign, of the non-native con- 
tacts relative to the native ones. 

We focus on 4 target conformations shown in Fig. 
1. There are 37 compact conformations (for the two- 
dimensional 16-monomer chain one has 69 compact 
4x4 conformations but for the Go-like model only 38 
of them remain different due to the end-to-end rever- 
sal symmetry and one is not accessible kinetically) 
which may act as the native conformations. Among 
these 5*1 and 5*4 shown in Fig. 1 are the two fastest 
folders at T ^ 0, whereas 5*2 and 6*3, also shown in 
Fig. 1, have intermediate folding properties. 
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FIG. 1. Four sequences studied in this paper . S\ and 
S3 have SE = 0. 

It should be noted that the repulsive non-native 
contact interaction are expected to improve the 



1 



foldicity ||^-|^ because they restrict the size of the 
relevant phase space. In our case, however, this ef- 
fect becomes so dramatic that for a > 0.1 sequences 
S'l and 53 can fold with a finite tfoid even at T = 0. 
Such an exotic phenomenon has been also observed 
in the HP model for some 13-monomer chains 0| 
so it is not restricted to the Go-like sequences that 
we study. Here, we study how this arises as a func- 
tion of a and show how do folding characteristics, 
i.e. characteristic temperatures and folding times, 
depend on a. The ability of S'l and 5*3 to fold at 
T = may be partly understood by the fact that for 
these sequences the repulsive non-native interactions 
reduce the number of local minima by two orders of 
magnitude compared to the case of a < 0. Fur- 
thermore, repulsive interactions dramatically affect 
partitioning of the phase space into regions associ- 
ated with the local energy minima. We demonstrate 
this by using the disconnectivity graph technique |^ 
and show, in particular, that connectivities to the 
folding funnel become simplified significantly. 
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FIG. 2. Temperature dependence of t/oid for S'l and 
52 and for 3 selected values of a, a = —0.4, and 1. The 
results are based on 2 - 4 batches, each containing 200 
trajectories from initial random conformations. 



|10[ . We obtain it, however, not only by the numeri- 
cal calculations for the lattice model but also by the 
analytical argument. The minimum folding time, 
imin, defined at the temperature where the folding 
is fastest has found to decrease with a and it gets 
saturated for a ^ oa. 

The folding dynamics of a chain is studied by a 
Monte Carlo procedure that satisfies the detailed 
balance condition |ll| , and was motivated by the 
studies presented in Ref . |0,|l2| . The dynamics allows 
for single and two-monomer (crankshaft) moves. For 
each conformation of the chain one has A possi- 
ble moves and the maximum value of A, Amax^ is 
equal to Amax = N + 2. In our 16-monomer case 
Amax = 18. For a conformation with A possible 
moves, probability to attempt any move is taken to 
be A/ Amax and probability not to do any attempt 
is 1-A/Amax In addition, probability to do a 

single move is reduced by the factor of 0.2 and to do 
the double move - by 0.8 |l^,@|. The attempts are 
rejected or accepted as in the standard Metropolis 
method. The folding time is equal to the total num- 
ber of Monte Carlo attemps divided by Amax ■ 

We have carried out the Monte Carlo simulations 
to determine the dependence of the median folding 
time, tfoid, on T and a. The results for Si and ^2 
are shown in Figure 2. For each temperature, tfoid 
is obtained based on 200 independent runs starting 
from random configurations. The results are aver- 
aged over 2-4 batches, of 200 trajectories each. 
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We have found that the folding temperature, Tf, 
increases with a. This result agrees with that of Ca- 
macho M for an effectively zero-dimensional model 



FIG. 3. Dependence of the folded fraction on a at 
T — for 4 sequences shown in Figure 1. The results 
are averaged over 4-6 simulations, each corresponding 
to 200 trajectories. 
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For sequence 5*2 the standard U-shape |]T3| for the 
temperature dependence of t fold is observed for all 
values of a. In other words, no qualitative change 
occurs if the non-native contact energies change from 
attractive to repulsive. In the case of sequence Si, 
however, for a = \ the standard U-shape disappears 
suggesting that the glass transition temperature Tg 
which is operationally defined as the value of the 
temperature at which the median time is equal to 
some cut-off value (usually this cut-off value is cho- 
sen to be 300000 Monte Carlo steps for the two- 
dimensional 16-monomer chain p^) becomes zero. 
In the standard scenario, at low temperatures the 
system may get trapped in some local minima and 
the folding process becomes extremly slow. The 
folding time is then governed by the Arrhenius law, 
tfoid ~ exp(5£'/T) [ p^|Jll|] , where 5E is the energy 
barrier energy (at T = one has tfoid oo). Thus, 
the absence of the fj-shape dependence suggests that 
the energy barrier 5E — 0. 
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FIG. 4. Dependence of T/, Tmin and tmin for sequence 
Si on a. The results are averaged over 4 batches, each 
containing 200 trajectories. 

In order to know whether 5E is exactly zero or 
not one has to study the folding at T = 0. At zero 
temperature the system gets trapped in some local 
minimum or in the native state. If the fraction of 
the trajectories from random conformations fold into 



the native state is bigger than 50%, then the chain 
is said to be folded and 5E = 0. 

The fraction of folded trajectories is shown in Fig. 
3 for S'l, 52, Sz and S/^. For Si and S^ this fraction 
becomes bigger than 50% for a > ac, where Ofc ~ 
0.1. Sequences 5*2 and 5*4 have SE for any value 
of a. It is interesting to mention that S4 folds even 
faster than S'l at T ^ but its foldicity becomes 
much worse at T = 0. Furthermore, the folding rate 
of 5*3 at T 7^ is slower than for the three other 
sequences and yet 5*3 can fold at T = 0. Thus the 
geometry of the native targets has a dramatic effect 
on the folding at T = 0. Among the 37 maximally 
compact 4x4 structures it is only S'l and S3 that 
do not obey the Arrhenius law. 
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FIG. 5. The same as in Figure 4 but for sequence ^2. 

We study the dependence of Tf on a for two typ- 
ical sequences Si and S2 shown in Fig. 1. The total 
number of conformations of the 16-monomer chain 
is only 802075 ||l^,|l^ and is amenable to exact enu- 
meration. This allows for an exact evaluation of the 
equilibrium parameters such as the folding temper- 
ature Tf. The latter is defined as a temperature 
at which the probability of occupancy of the native 
state is 1/2. The results for sequences Si and S2 are 
shown in Fig. 4 and 5. Obviously, Tf increases with 
a but this dependence gets weaker for larger values 
of a. 
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The increase of Tf with a may be understood in 
the foUowing simple way. Let Ptq (a) be the proba- 
bihty of occupying the native state of conformation 
To. Then 



Pro (a) 



exp' 



(2) 



where Eq is the energy of the native state, n is the 
number of non-native contacts and is the part of 
energy in conformation T which corresponds to the 
native contacts. Then 



da 



> (3 <n>T>0 



(3) 



where < n >t is the average number of non-native 
bonds at temperature T. So the probabiUty of be- 
ing in the native state cannot decrease with a. Tj, 
therefore, should increase with increasing a and then 
become a-independent. 




We now focus on the a dependence of the tem- 
perature at which the folding time is minimal, Tmm- 
The results for and 52 are shown in Fig. 4 and 
Fig. 5. For positive values of a, Tf and Tmin sue 
comparable for both sequences and they should be 
good folders For large negative values of a, 

Tmin is bigger than Tf but 5*1 and 52 remain good 
folders by the Thirumalai-Camacho criterion [ [l7||lq | 
- the peaks of the structural susceptibility and the 
specific heat coincide for this interval of a. (For 
a = 0, the coincidence of the peaks is, in fact, a gen- 
eral feature of the Go models because the proximity 
to the native state and occurence of rapid changes 
in an average energy as a function of T are both 
controlled by establishment of the same native con- 
tacts). 
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FIG. 7. The disconnectivity graph for Si and a — 0. 
Nc is a symbolic notation for the index labelling the local 
energy minima. 



FIG. 6. Dependence of Nim on a for Si and 5*2. 

It should be noted that in Camacho's model ||] 
Tf was found to increase with a {a < 0) linearly. 
Our results presented in Figures 4 and 5 show that 
the region of a where one can observe the linear 
dependence is rather narrow. Such region becomes 
much wider, for example, in the case of the 27- 
monomer chain in three dimensions (the results are 
not shown). Overall, the results for Tf, shown in Fig. 
4 and 5, demonstrate that the non- native contact 
repulsive energies improve both the thermodynamic 
stability and dynamical characteristics of folding. 



Fig. 4 and 5 also show the dependence of the mini- 
mal folding time, tmim on a. The dependence seems 
to saturate at large values of a. If one extends the 
Camacho result to positive a then tmin should de- 
crease with a exponentially Our results suggest 
that such conclusion is valid only for a < but not 
for a > 0. Thus, the repulsive non-native contact 
energies make the polypeptide chain to fold faster 
which is similar to what has been observed in Ref. 

i|. 

In order to understand why the repulsion improves 
the folding so much we study the dependence of 
number of local minima, iV/,„, on a. The results for 
Si and 52 are shown in Fig. 6. Nim of 5i is found to 
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be smaller than for 82- Clearly, the number of local 
minima strongly depends on a and for a > it be- 
come by 2 orders of magnitude smaller than that for 
a < 0. So, the positive values of a make the energy 
landscape less rugged and the folding dynamics get 
faster. 



the local energy minima that have links to the na- 
tive state below any predetermined energy thresh- 
old. Thus, the comparison of the disconnectivity 
graphs also shows that the repulsion may facilitate 
the folding substantially. 
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FIG. 8. The same as in Fig. 7 but for a = 1. 

In order to get more insight into the nature of the 
energy landscapes in the models studied here, we use 
the disconnectivity graph technique which maps the 
potential energy surface onto the set of local minima 
[^,|9|. The technique involves checking what local 
energy minima are connected by pathways (sets of 
moves that are allowed kinetically) that do not ex- 
ceed a given total threshold energy. For each value 
of this energy the minima are divided into discon- 
nected sets of mutually accessible minima separated 
by barriers. The local minima which share the low- 
est energy threshold are joined at a node by lines 
and are called a basin corresponding to the thresh- 
old. The procedure of construction is stoped when 
one gets only one basin for all of the minima. For Si 
the number of local minima is equal to 152 and 81 for 
a — Q and a — \, respectively. The construction of 
the graphs, therefore, may be done exactly. The dis- 
connectivity graphs obtained for Si and a = and 
a — 1 are shown in Fig. 7 and Fig. 8, respectively. 
For both values of a sequence Si is a good folder 
and consequently, the structure corresponding to a 
folding funnel is clearly visible. However, the funnel 
for the repulsive case of a = 1 has a significantly less 
complex pattern than for a ~ 0. Moreover, the re- 
pulsive non-native interactions reduce the number of 
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FIG. 9. The density of the (V-shape) local minima, 
Plm, in which the system gets trapped at T = is 
plotted versus the number of native contacts for 5*1 and 
5*2. We choose a — 0.1 and a = —0.1. The results 
are obtained for the batch of 200 trajectories. In the 
case of the local minima are different for a = —0.1 
and Q = 0.1. For S2 there are 19 common local minima 
which are marked by dotted lines. 

We now address the question of what happens 
with trapped local minima when a changes from 
negative to positive values. Fig. 9 shows the his- 
togram of the local minima for a = —0.1 and 0.1 
for sequence 5*1 and 82- Interestingly, for 5*1 none 
of the local minima obtained for a = —0.1 appears 
for a = 0.1. The situation changes dramatically for 
5*2 for which 19 local minima are common for both 
a — —0.1 and a = 0.1. Our results suggest that for 
the sequence with SE — the local minima in which 
the chain is trapped at the negative values of a are 
effectively avoided if a is changed to positive values. 

In conclusion, we state that for the simple ex- 
tended Go model Tf (tmm) increases (decreases) 
with the non-native contact energy and it gets sat- 
urated for large values of a. The complexity of the 
disconnectivity trees becomes reduced on making a 
more and more repulsive so that some two dimen- 
sional sequences may even lose the Arrhenius like 
behavior of tfoid at low temperatures. It would 
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be interesting to determine whether there are any 

three-dimensional Go-Ukc sequences that fold even 
at T = 0. Another important question is what kind 
of an effective parameter a, or its conceptual equiv- 
alent, characterizes real proteins. 
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